Thermodynamic Formalism of the Harmonic Measure of Diffusion Limited 
Aggregates: Phase Transition and Converged /(a) 
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We study the nature of the phase transition in the multifractal formalism of the harmonic measure 
of Diffusion Limited Aggregates (DLA). Contrary to previous work that relied on random walk 
simulations or ad-hoc models to estimate the low probability events of deep fjord penetration, we 
employ the method of iterated conformal maps to obtain an accurate computation of the probability 
of the rarest events. We resolve probabilities as small as 10" . We show that the generalized 
dimensions D q are infinite for q < q* , where q* — — 0.17± 0.02. In the language of f(a) this means 
that a max is finite. We present a converged f(a) curve. 
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Since its introduction in 1981 W the model of Diffu- 
sion Limited Aggregation (DLA) has posed a challenge to 
our understanding of fractal and multifractal phenomena. 
DLA is a paradigmatic example for the spontaneous gen- 
eration of fractal objects by simple dynamical rules (be- 
ing generated by random walkers); its harmonic measure, 
which is the probability for a random walker to hit the 
surface, had been one of the first studied examples of mul- 
tifractal measures outside the realm of ergodic measures 
in dynamical systems B. The multifractal properties 
stem from the extreme contrast between the probability 
to hit the tips of the DLA compared with penetrating 
the fjords. 

The multifractal properties of the harmonic measure 
of the DLA are conveniently studied in the context of 
the generalized dimensions D q , and the associated f{a) 
function ||Q . The simplest definition of the generalized 
dimensions is in terms of a uniform covering of the bound- 
ary of a DLA cluster with boxes of size £, and measuring 
the probability for a random walker coming from infinity 
to hit a piece of boundary which belongs to the i'th box. 
Denoting this probability by Pi(£), one considers M 



D„ = lim 



i log £< *?(*) 



»o q — 1 log £ 



(1) 



It is well known by now that the existence of an interest- 
ing spectrum of values D q is related to the probabilities 
Pi(t) having a spectrum of "singularities" in the sense 
that Pi(£) ~ £ a with a taking on values from a range 
amin < ol < a max . The frequency of observation of a 
particular value of a is determined by the the function 
f(a) where (with r(q) = (q — l)D q ) 



/(a) = aq(a) -rfg(a) 



dr(q) 



a(q) 



(2) 



The understanding of the multifractal properties and 
the associated f(a) spectrum of DLA clusters have been 
a long standing issue. Of particular interest are the val- 
ues of the minimal and maximal values, a m i n and a max , 



relating to the largest and smallest growth probabilities, 
respectively. As a DLA cluster grows the large branches 
screen the deep fjords more and more and the probability 
for a random walker to get into these fjords (say around 
the seed of the cluster) becomes smaller and smaller. A 
small growth probability corresponds to a large value 
of a. Previous literature hardly agrees about the ac- 
tual value of a max . Ensemble averages of the harmonic 
measure of DLA clusters indicated a rather large value 
of a max ~ 8 pi. In subsequent experiments on non- 
Newtonian fluids P and on viscous fingers J7]], similar 
large values of a max were also observed. These numeri- 
cal and experimental indications of a very large value of 
otmax led to a conjecture that, in the limit of a large, self- 
similar cluster some fjords will be exponentially screened 
and thus causing a max — ► oo H. 

If indeed a max — > oo, this can be interpreted as a phase 
transition |9J (non-analyticity) in the q dependence of D q , 
at a value of q satisfying q > 0. If the transition takes 
place for a value q < then a max is finite. Lee and Stan- 
ley pi proposed that a max diverges like R 2 /lnR with 
R being the radius of the cluster. Schwarzer et al. |ll]] 
proposed that a max diverges only logarithmically in the 
number of added particles. Blumenfeld and Aharony |1J] 
proposed that channel-shaped fjords are important and 
proposed that a max ~ f ^ , f w here M is the mass of the 
cluster; Harris and Cohen |l3|, on the other hand, argued 
that straight channels might be so rare that they do not 
make a noticeable contribution, and a max is finite, in 
agreement with Ball and Blumenfeld who proposed pj] 
that a max is bounded below 11. Obviously, the issue was 
not quite settled. The difficulty is that it is very hard to 
estimate the smallest growth probabilities using models 
or direct numerical simulations. 

In this Letter we use the method of iterated conformal 
maps to offer an accurate determination of the proba- 
bility for the rarest events. We propose that using this 
method we can settle the issue in a conclusive way. Our 
result is that a max exists and the phase transition occurs 
at a q value that is slightly negative. In this method one 



studies DLA by constructing 3>( n \w) which conformally 
maps the exterior of the unit circle e in the mathe- 
matical w-plane onto the complement of the (simply- 
connected) cluster of n particles in the physical z-plane 
Eq O . The unit circle is mapped onto the boundary of 
the cluster. The map <&( n >(w) is made from compositions 
of elementary maps </>A,e, 



*(»)(«,) = $("-!) (cf> Xn , 6n (w)) 



(3) 



where the elementary map (f>\ t g transforms the unit cir- 
cle to a circle with a semi-circular "bump" of linear size 
vA around the point w = e l8 . We use below the same 
map </>a,0 that was employed in p"5| |i~9|] . With this map 
<£>(«) ( w ) adds on a semi-circular new bump to the im- 
age of the unit circle under ^ n ~ 1 \w). The bumps in 
the z-plane simulate the accreted particles in the phys- 
ical space formulation of the growth process. Since we 
want to have fixed size bumps in the physical space, say 
of fixed area An , we choose in the nth step 



A„ = 



An 



|$(n-l)'( e i0„)|5 



(4) 



The recursive dynamics can be represented as iterations 
of the map 4>\ n ,e n {w), 



$ (n) (w) = (f>\ u e 1 o <t>\ 2 ,e 2 o . . . o 4> Xnfin (uj) 



(5) 



It had been demonstrated before that this method repre- 
sents DLA accurately, providing many analytic insights 
that are not available otherwise [T^Jl9|. For our purposes 
here we quote a result established in |16(| , which is 

(XI) = (1/2tt) r K(0)d6 ~ n-^ D ^' D . (6) 



To compute r(q) we rewrite this average as 



(K) 



ds 



dd 
ds 



K(s) 



X q+1/2 (s) 
'Is - J S) , (7) 



fM 



where s is the arc-length of the physical boundary of 
the cluster. In th e last equality we used the fact that 
\d9/ds\ = y/X n /X . We stress at this point that in order 
to measure these moments for q < we must go into 
arc-length representation. 

To make this crucial point clear we discuss briefly what 
happens if one attempts to compute the moments from 
the definition (|6|). Having at hand the conformal map 
<j>(")(e l6 '), one can choose randomly as many points on 
the unit circle [0, 2ir] as one wishes, obtain as many (ac- 
curate) values of A n , and try to compute the integral as 
a finite sum. The problem is of course that using such 
an approach the fjords are not resolved. To see this we 
show in Fig.l panel a the region of a typical cluster of 50 
000 particles that is being visited by a random search on 
the unit circle. Like in direct simulations using random 



walks, the rarest events are not probed, and no serious 
conclusion regarding the phase transition is possible. 
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FIG. 1. panel a: the boundary of the cluster probed by a 
random search with respect to the harmonic measure. Panel 
b: the boundary of the cluster probed by the present method. 

Another method that cannot work is to try to com- 
pute by sampling on the arc-length in a naive way. The 
reason is that the inverse map [<&(")] -1 (s) cannot resolve 
8 values that belong to deep fjords. As the growth pro- 
ceeds, reparametrization squeezes the 9 values that map 
to fjords into minute intervals, below the computer nu- 
merical resolution. To compute the values of X n (s) ef- 
fectively we must use the full power of our iterated con- 
formal dynamics, carrying the history with us, to iterate 
forward and backward at will to resolve accurately the 9, 
X values of any given particle on the fully grown cluster. 

To do this we recognize that every time we grow a semi- 
circular bump we generate two new branch-cuts in the 
map <I > ( n ' . We find the position on the boundary between 
every two branch-cuts, and there compute the value of 
A n . The first step in our algorithm is to generate the 
location of these points intermediate to the branch-cuts 
pTj. Each branch-cut has a preimage on the unit circle 



which will be indexed with 3 indices, 



k{l) _ 
°J4 = 



expp 



Mth 



The index j represents the generation when the branch- 
cut was created (i.e. when the jth particle was grown). 
The index £ stands for the generation at which the anal- 
ysis is being done (i.e. when the cluster has £ particles). 
The index k represents the position of the branch-cut 
along the arc-length, and it is a function of the gen- 
eration £. Note that since bumps may overlap during 
growth, branch-cuts are then covered, and therefore the 
maximal k, k max < 2£. After each iteration the preimage 
of each branch-cut moves on the unit circle, but its phys- 
ical position remains. This leads to the equation that 
relates the indices of a still exposed branch-cut that was 
created at generation j to a later generation n: 
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-l^ 3,3 



>) 



(8) 



Note that the sorting indices k(j) are not simply related 
to k(n), and needs to be tracked as follows. Suppose 
that the list w, T "J7i is available. In the nth genera- 
tion we choose randomly a new 9 n , and find two new 
branch-cuts which on the unit circle are at angles 0„ . If 
one (or very rarely more) branch-cut of the updated list 

cf>^ e {Wjn-i ) i s covered, it is eliminated from the list, 
and together with the sorted new pair we make the list 
w i n ■ Having a cluster of n particles we now consider all 

neighboring pairs of preimages w, ™ and Wj™ , that 
very well may have been created at two different gener- 
ations j and J. The larger of these indices (J without 
loss of generality) determines the generation of the inter- 
mediate position at which we want to compute the field. 
We want to find the preimage Uj y of this mid-point on 
the unit circle , to compute Xkln) there accurately. Using 
definition (JS|) we find the preimage 

arg(ug>) = [arg(«,jy>) + arg(<^ +1 )]/2 . (9) 

In Fig.l panel b we show, for the same cluster of 50 000, 
the map ^'"("jj ) with k(n) running between 1 and 
fcmax, with J being the corresponding generation of cre- 
ation of the mid point. We see that now all the particles 
are probed, and every single value of Xk( n ) can be com- 
puted. However, to compute these Xk( n ) accurately, we 
define (in analogy to Eq.(0)) for every J < rn < n 



k(r. 



J,m — n„,«„ ° • ' • ° ( f\ J + u 9 J+1 \ u J,J ) 



(10) 



Finally X k t n \ is computed from the definition (H) with 

x ^'(u*}?) . (11) 

This calculation is optimally accurate since we avoid as 
much as possible the effects of the exponential shrink- 
ing of low probability regions on the unit circle. Each 



derivative in (nil) is computed using information from a 
generation in which points on the unit circle are optimally 
resolved. 
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FIG. 2. Panel a: the calculated function r(q) for clusters of 
n particles, with n =10 000,15 000, 25 000 and 30 000. Panel 
b: the second derivative of r(q) with respect to q. 



The integral (Q) is then estimated as the finite sum 
\^oEi(n)^L)- We should stress that for clusters of 
the order of 30 000 particles we already compute, us- 
ing this algorithm, Aj,(„) values of the order of 10~ 70 . 
To find the equivalent small probabilities using random 
walks would require about 10 attempts to see them just 
once. This is of course impossible, explaining the lasting 
confusion about the issue of the phase transition in this 
problem. This also means that all the f(a) curves that 
were computed before [plplj did not converge. Note that 
in our calculation the small values of Afc(„j are obtained 
from multiplications rather than additions, and therefore 
can be trusted. 

Having the accurate values Xki n -j we can now compute 



the moments (|6|). Since the scaling form on the RHS 
includes unknown coefficients, we compute the values of 
r(q) by dividing (A&) by (A|), estimating 



r(q) » -D 



log(A«)-log(A|) 
log n — log n 



(12) 
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Results for r(q) for increasing values of n and n are shown 
in Fig. 2, panel (a). It is seen that the value of r(q) 
appears to grow without bound for q negative. The exis- 
tence of a phase transition is however best indicated by 
measuring the derivatives of r(q) with respect to q. In 
Fig. 2 panel b we show the second derivative, indicating 
a phase transition at a value of q that recedes away from 
q = when n increases. Due to the great accuracy of our 
measurement of A we can estimate already with clusters 
as small as 20-30 000 the q value of the phase transition 
to q = — 0.17±0.02. The fact that this value is very close 
to the converged value can be seen from the f(a) curve 
which is plotted in Fig. 3. A test of convergence is that 
the slope of this function where it becomes essentially lin- 
ear must agree with the q value of the phase transition. 
The straight line shown in Fig. 3 has the slope of -0.17, 
and it indeed approximates very accurately the slope of 
the f(a) curve where it ends and stops being analytic. 
The reader should also note that the peak of the curve 
agrees with D rs 1.71, as well as the fact that r(3) is also 
D as expected in this problem. The value of a max is close 
to 20, which is higher than anything predicted before. It 
is nevertheless finite. We believe that this function is well 
converged, in contradistinction with past calculations. 




FIG. 3. The calculated function f(a) using r(q) calculated 
from a cluster with n = 30 000 particles. This f(a) is almost 
indistinguishable from the one computed with n = 25000 par- 
ticles. We propose that this function is well converged. The 
black dot denotes where the curve ends, being tangent to the 
line with slope -0.17. 
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